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ABSTRACT 

Kepler-78b is a transiting planet that is 1.2 times the radius of Earth and orbits a young, active K 
dwarf every 8 hours. The mass of Kepler-78b has been independently reported by two teams based on 
radial velocity measurements using the HIRES and HARPS-N spectrographs. Due to the active nature 
of the host star, a stellar activity model is required to distinguish and isolate the planetary signal in 
radial velocity data. Whereas previous studies tested parametric stellar activity models, we modeled 
this system using nonparametric Gaussian process (GP) regression. We produced a GP regression of 
relevant Kepler photometry. We then use the posterior parameter distribution for our photometric 
fit as a prior for our simultaneous GP + Keplerian orbit models of the radial velocity datasets. We 
tested three simple kernel functions for our GP regressions. Based on a Bayesian likelihood analysis, we 
selected a quasi-periodic kernel model with GP hyperparameters coupled between the two RV datasets, 
giving a Doppler amplitude of 1.86 ± 0.25 m s“^and supporting our belief that the correlated noise 
we are modeling is astrophysical. The corresponding mass of l-87lo!26 ^0 consistent with that 
measured in previous studies, and more robust due to our nonparametric signal estimation. Based 
on our mass and the radius measurement from transit photometry, Kepler-78b has a bulk density 
of 6.01};4 g cm“^. We estimate that Kepler-78b is 32 ± 26% iron using a two-component rock-iron 
model. This is consistent with an Earth-like composition, with uncertainty spanning Moon-like to 
Mercury-like compositions. 


1. INTRODUGTION 


Measuring the radial velocity (RV) of a planet’s host 
star is the most common method to measure the mass 
of a planet. The planetary RV signal, or Doppler ampli¬ 
tude, is directly related to the planet mass. If we know 
the mass of the star, we can find the mass of the planet 
from its Doppler amplitude. If the planet radius is also 
known, we can then calculate the planet’s density and es¬ 
timate a composition. Precise density and composition 
information is available only for a handful of transiting 
rocky planets. 


However, stellar activity can produce spurious RV sig¬ 
nals larger than some planetary signals. In young, spot¬ 
ted stars, this activity can cause RV variation s much 


larger than most observed planetary RV signals (Hillen- _ 
brand et al.|2Q15 ), and even in relatively quiet stars such 


as the bun, these variations are at least one order of 


magnitude larger than the RV signal of Earth (Meunier 


et al. 2010). Before we can have any hope of confirming 


the discovery of an Earthlike planet around any star, we 
must account for the contribution of stellar activity. 


Discovered in 2013, the ra dius of Kepler-78b was mea - 
sured to be 1.16 ± 0.16 R 0 (|Sanchis-Ojeda et al.||2013|). 
Subsequently, the mass of the planet was measured and 
reported simultaneously by two teams. Based on the ra- 
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bution from correlated noise activity while measuring the 
planetary Doppler signal. H13 used a sum of sinusoids 
at the stellar rotation period and its harmonics to model 
the correlated noise activity, accounting for rotational 
variability while neglecting any component of the pre¬ 
dominant signal that was not periodic at precisely the 
rotation period of the star or its aliases. 

P13 noted the fact that the planetary and stellar RV 
signal timescales differ by over an order of magnitude: 
the orbital period of the planet is only 8.5 hours, whereas 
the stellar activity-induced signals are modulated by 
the stellar rotation period of 12.8 days. This allowed 
P13 to test a floating chunk offset model in which the 
free parameters were the planetary Doppler amplitude 
and an RV zero point to represent the noise signal, as¬ 
sumed constant for each night of observations. How¬ 
ever, P13 found that the evidence ratio for their float¬ 
ing chunk offset model was significantly lower than the 
sum of sinusoids model, given the large number of free 
parameters, and used the parametric, si nusoidal model 
for their reported mass measurement. [Howard et al. 
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(2013)used Keck/HIRES and measur ed the mass of Ke¬ 
pler 78b to be 1.69 ± O.41M0, while Pepe et al. (2013) 
used TNG/HARPS-N and measured a mass of 1.86^0.25 
M 0 . Despite the difference in their methods and obser¬ 
vational techniques, the H13 and P13 planet mass mea¬ 
surements are consistent, attesting to the robustness of 
the results. 


Further exploration into nonparametric estimation is 
justified by the limited scope of the purely parametric 
models tested previously. We draw attention to the com¬ 
parable case of HST transmis sion spectroscopy o f hot 
Jupiter HD 1897 33, s tudied by Swain et al. 


son 


^ al. (2011) and lGibson et al? |(|2U12|). |Gibson et al. 
^ note that botn previous studies of the transmis- 


(| 2 U 12 p note tiiat botii previous 
Sion spectroscopy used linear basis functions to account 
for systematic errors, and argue that this is not sufficient 
to account for instrumental systematics, and therefore 
provides an unrealis tic treatment of the uncertainties. 


Gibson et al. ( 2012 ) reanalyzed the spectroscopic data 
with a Gaussian process to marginalize ignorance of the 
functional form of the systematics, giving larger but more 
robust errors. Similarly, we have reanalyzed the RV data 
of Kepler-78b with several Gaussian process estimators 
to find a model of the RVs with the highest evidence, 
resulting in more robust uncertainties of the mass, den¬ 
sity, and composition of Kepler-78b. In this case, we are 
using Gaussian process to describe a signal of astrophys- 
ical origin, rather than instrumental, which we confirm 
through tests described in our Results section. 


In this study, we combine the RV measurements from 
H13 and P13. We describe each radial velocity dataset 
with a Gaussian process (GP) regression combined with 
a Keplerian orbit signal at the orbital period and phase 
of the planet (both known precisely from the transit pho¬ 
tometry). GP regression is a nonparametric method for 
modeling correlated RV noise, and can provide a more ro¬ 
bust noise model because of its flexibility over parametric 
models. The observations analyzed are described in §2. 
Our analysis of the combined RV dataset is described in 
depth in §3, starting with a review of Gaussian process 
regression and the benefit of a nonparametric model in 
§ 3.1-3.3, and a focus on how different GP models were 
tested and chosen for the data in sections § 3.4 and 3.5. 
In §4, we report our results for the planet Doppler am¬ 
plitude and mass, and discuss the selection of the cho¬ 
sen analysis model. Using the known properties of the 
host star (summarized in Tableand the planet radius, 
we recalculate the planetary density. In §5, we compare 
our mass measurement and density calculation to pre¬ 
vious results and discuss possible composition scenarios 
for Kepler-78b. We explore the applications of this tech¬ 
nique to other RV datasets in § 6 . 


2 . OBSERVATIONS 

In order to determine the Doppler amplitude of Kepler- 
78b, it was necessary to first confirm the general struc¬ 
ture of the stellar activity by modeling the Kepler pho¬ 
tometry. These model parameters were then used to pro¬ 
vide reasonable initial hyperparameters for the RV anal¬ 
ysis. The observations taken are described below. 


TABLE 1 

Stellar Properties 


Property 

Value 

Name 

Kepler-78/KIC 8435766/Tyclio 3147-188-1 

Age 

625 ± 150 million years 

Ysin{i) 

2.6 ± 0.5 km s“^ 

Mass, Mstar 

0.83 ± 0.05 M© 

Inclination, i 

75 . 2 II® degrees 

Rotation period Prot 

12.5 days 


Note. — The values in Table [T] are have been taken from lHoward eti 
[£T1 ( [2QT^ . 


2 . 1 . Kepler photometry 

The Kepler telescope obtained photometry of approx¬ 
imately 150,000 objects in the Kepler field for its four 
year lifetime from 2009 to 2013. Photometry of Kepler- 
78 was gathered at a 30-minute cadence. In this study, 
we train a GP on the photometric light curve of Kepler- 
78 to determine the evolution and rotation timescales of 
the stellar activity, hyperparameters of the GP regres¬ 
sion kernel. However, as shown in Equation a matrix 
inversion is required to calculate the log posterior likeli¬ 
hood of a GP kernel with a given set of hyperparameters, 
a computationally intensive process with compute time 
proportional to N^. Therefore, it was necessary to rebin 
the Kepler photometry to one point every 5 hours (aver¬ 
aging every ten points together) in order to find the best- 
fit hyper parameters of a full quarter of photometry with 
our computing resources. While this does marginalize 
over the planetary signal, we are only using the photom¬ 
etry to estimate the underlying components of the stel¬ 
lar activity which vary on timescales much longer than 5 
hours, and the photometric effect of the planet is small 
compared to that of the stellar activity. This photome¬ 
try, taken from the 16th quarter of Kepler observations, 
is not concurrent with the RV measurements; however, 
since it was taken only ^55 days before the RV mea¬ 
surements, we can assume that the variation seen in the 
photometry is at least related to the RV variation, and 
would exhibit similar structure in the temporal dimen¬ 
sion. Figure illustrates the binned photometric mea¬ 
surements as well as the GP regression with the best-fit 
kernel hyperparameters to the photometry. 

2 . 2 . HIRES and HARPS-N RV data 

The RV observations were obtained shortly after the 
malfunctioning of KeplerA reaction wheels in summer 
2013, which prevented the high level of photometric pre¬ 
cision achieved in the first era of the spacecraft’s life¬ 
time and required that the telescope point at new tar¬ 
gets. H13 observed Kepler-78 with HIRES on the Keck 
I telescope on Manna Kea. They obtained eight nights 
of data during each of which approximately 9 or more 
measurements were obtained: three measurements were 
taken at successive thirty minute intervals at three sepa¬ 
rate occasions per night. In addition, five individual mea¬ 
surements were taken on separate nights after the rest of 
the data were collected to provide a longer baseline over 
which to evaluate the stellar activity. The dataset con¬ 
sists of 84 RVs over 45 nights. The P13 team observed 
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Kepler-78 with HARPS-N on the Telescopio Nazionale 
Galilei in the Canary Islands. They were able to obtain 
109 measurements over 97 nights. Most of their obser¬ 
vations were made in a consecutive six-day period. On 
subsequent nights only one or two measurements were 
made per night, giving a total time baseline of approx¬ 
imately three months. Figure illustrates the HIRES 
and HARPS-N measurements with measurement errors 
shown, with GP regressions overplotted. 


3. METHODS 


We model the photometric variations observed by Ke¬ 
pler with a GP regression with a quasi-periodic kernel. 
We then model the radial velocity measurements with a 
quasi-periodic kernel GP regression, using the kernel pe¬ 
riod hyperparameter, 6>phot: fo train the RV GP period 
hyperparameter 6 >rv. This assumes that the predomi¬ 
nant signal observed in both datasets comes from the 
star, and thus will modulate at the stellar rotation pe¬ 
riod. When we model the photometric data, we recover a 
period hyper parameter that agrees with o ther autocorre¬ 
lation anal yses of the Kepler photometry (Sanchis-Ojeda 
et al.|2Q13 ), supporting this assumption. We create a GP 
regression of each RV dataset with a quasi-periodic ker¬ 
nel, where the period hyperparameter has been trained 
on the period hyperparameter of the photometric regres¬ 
sion. We then fit all the RV GP hyperparameters, h, 
6>, re, and A in the quasi-periodic case, a, a white noise 
parameter, as well as the Doppler amplitude K of a. Ke- 
plerian signal at the known orbital period and phase of 
the planet, to the RV data simultaneously. We find that 
all of our GP + Keplerian models are able to recover 
the planetary Doppler amplitude K at a value consistent 
with previous work. 


3.1. Gaussian Proeess Regression: Coneept 


spot models sueh as those used by HI3 and P13 gave in- 
eorreet masses and uneertainties. Thus, it is important 
to test many time series teehniques, and further explore 
the novel applieation of GPs in Doppler analys is. We ex¬ 
pand upon t he GP training technique used in |Haywood| 


et al 


(2014) here. 


3.2. Gaussian Proeess Govarianee Kernel Ghoiee 

Einding the best GP regression requires choosing a 
kernel and initial hyperparameters, evaluating the likeli¬ 
hood of those hyper parameter values, and then iterating 
through parameter space until the most likely values are 
found. The squared exponential kernel, for example, de¬ 
fines a covariance matrix through an operator. 


T^ij = k{ti,tj) = h^exp 


(^) ■ 


( 1 ) 


where h is the covariance amplitude, and A the covariance 
length scale. The amplitude observed is described by h, 
while A is a characteristic timescale over which the data 
is going to be correlated. 


We discuss other GP kernels and the inferred physical 
meaning of their hyperparameters in Table 

The logarithm of the posterior likelihood of the GP 
regression is calculated as 


log[£(r)] = - hog|S| - |log(27r), (2) 

where r is the vector of residuals after removal of the 
(optional) mean function, 5] is the covariance matrix, 
and n the number of data points. A prior term, >Cprior, 
can be added to the likelihood to account for any priors 
placed on the hyperparameters. Eor example, we apply 
the Gaussian prior 


Gaussian process is a nonparametric method to de¬ 
scribe a dataset by evaluating correlations between n 
data points through a covariance kernel. This kernel 
describes the relationship of each point in the dataset 
to each other point, and can be expressed as an n x n 
matrix (subsequently referred to as the covariance ma¬ 
trix). The kernel is a function of hyperparameters. More 
complicated kernels can have more hyperparameters that 
characterize different qualities of the correlations in the 
data, such as various periods, characteristic amplitudes 
and length scales, ete. 


Gaussian process re gression is widely used in the field 
of machine learning (Neal 1997| |Herbrich et al. 2003 


Quihonero-Gandela V Rasmussenj |2005[ |Wang et^ 

2008|. Gibson et al. (|2012 ) introduced the technique to 
the held of exoplanets through analysis of transmission 
spectroscopy to model correlated noise in the instrumen¬ 
tal systematics of HST/NI CMOS, as described i n §1. 
Goncurrently with this work, Haywood et al. (2014) have 
demonstrated the technique of GF modeling of RV and 
photometric signals for the GoRoT-7 planetary system, 
first modeling the photometry with a GP and then us¬ 
ing the photometric GP hyperp arameters to train the 
initial RV GP hyperparameters. \Haywood et aL\ l[2014\l 
demonstrated that in the ease of CoKoT-7b, parametrie 


_ 1 

C ■ =e ^ 

^prior ^ 

to restrict the hyperparameter 0 for the RV data, as the 
period of the correlated noise signal is equivalent to the 
stellar rotation period Prot in both the photometric and 
RV datasets. We restrict the period hyper parameter 0 
in the RV regressions within a Gaussian of width cfq de¬ 
termined by the posterior distribution of the period hy¬ 
perparameter found for our photometric GP regression. 
Prior knowledge of this hyperparameter helps to ensure 
convergence of the other hyperparameters of the RV re¬ 
gression, as the RV data is not as well sampled as the 
photometry. 




( 3 ) 


This likelihood calculation can be used to identify the 
best fit GP hyperparameters. Eor more details on how 
the covariance kernel and parameter boundary conditions 
were chosen, refer to § 4.1. Eor a more complete descrip¬ 
tion of Gaussian pr ocess regression and posterior likeli- 
hood evaluation, see Rasmussen & Williams (2006). 


3.3. Gaussian Proeess and Stellar Aetivity 

The presence of magnetic features on the stellar sur¬ 
face induces variation in a star’s RV that can mimic a 
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Fig. 1.— Flux of Kepler-78 versus time, as recorded by the Kepler spacecraft during its 16th quarter of observations. The data are plotted 
as black points, with photometric errors shown. The blue line corresponds to the Gaussian process regression of the photometry with the 
best-fit kernel hyperparameters, and the shaded regions correspond to l-cr and 2-a uncertainties, as defined by the posterior distributions 
of the kernel hyperparameters. The Kepler data has been shown as fitted with a single quasi-periodic kernel function and a white-noise 
stellar jitter parameter. 


TABLE 2 

Gaussian Process Kernel Options 


Name 


Squared exponential 

Periodic 


Mathematical expression 


h^exp 
h^exp 




sin^ —tj )/9] 
2^2 


Quasi-Periodic 


h^exp 


sin^ —tj )/9] 

2w‘^ 



Hyperparameters^ 

Gomments 

/i, A 

h amplitude of covariance function, 

A a characteristic timescale 

/l, 0, w 

0 equivalent to Prot, 

w represents coherence scale, similar to A expressed 
as a fraction of 9 dependent on recurrent features 

/i, 0 , r;, a 

w coherence scale tied to periodic variation 

while characteristic timescale A tied to aperiodic variation. 


Note. — The name of kernel functions and hyperparameters in Table are taken from [Rasmussen Williams| ( |2006| ). 

Each kernel X)ij can be modified to include an additional hyperparameter, a white noise term cr^ by adding one in quadrature: cr^li. 


planet’s orbital RV signal. Starspots are regions of high 
magnetic fields on the surface of a star that appear rel¬ 
atively cooler and darker than the surrounding photo¬ 
sphere. As the starspots move across the disk of a star 
(due to stellar rotation), the flux balance between the 
redshifted and blueshifted halves of the star is broken, 
producing a shift in the centroid of the spectral lines. 
This shift corresponds to an apparent RV signal and 


change as the spots move across the stellar disk (Du- 
musque et al. 2Q11[). Starspots, as well as networks of 


strongly magnetized flux tubes known as faculae, inhibit 
the convective processes taking place on the stellar sur¬ 
face, thus reducing the net blueshift produced by the up 
flow of hot, bright granules. This h as been shown to be 
the dominant RV effect on the Sun (Meunier et a l. 12010) 
and on other Sun-like stars (Haywood et al. |2U14[ ). These 
activity-induced RV signals are modulated by the stellar 
rotation and the surface features evolve with time, re¬ 
sulting in a quasi-periodic, RV signal. The lifetimes of 
surface features are poorly constrained, making this ef¬ 
fect difficult to model parametrically. The nonparamet- 
ric nature of the GP regression provides an ideal frame¬ 
work for studying stellar surface features in RV. We can 
choose a GP kernel that reflects the frequency structure 
of the stellar activity. The resultant GP model is flexible 
enough to account for the evolution of magnetic features 


on the stellar surface while keeping a statistical “mem¬ 
ory” of the stellar activity patterns. 


3.4. Photometric Model 

We fit the photometric data using GP regressions with 
three different kernels. We test a squared exponential, 
periodic and quasi-periodic kernel (Table [^. We test 
kernels with and without an additional white noise term 
added in quadrature to the likelihood (Eq. 2). The best- 
fit kernel hyper parameters are found via Markov Chain 
Monte Carlo analysis powered by t he emcee Python 
package (Foreman-Mackey et al. 2013), described in more 
detail in the following section. We select the quasi- 
periodic model with a white noise term based on our as¬ 
sumption that the predominant signal in both datasets is 
tied to the quasi-periodic variation of the stellar surface, 
and confirm that the quasi-periodic model is appropriate 
through a visual check and a reduced analysis. We 
have plotted the data and GP regression with the best- 
fit kernel hyperparameters in Figure We ensure that 
these hyperparameters are robust by performing a GP 
regression to the photometric data after less stringent 
binning, and recovering the same best-fit kernel hyper¬ 
parameter values within errors. 






























Kepler-78b Mass Measurement with Gaussian Process 


5 


3.5. RV Model 


We fit our photometric and RV datasets with the same 
kernel operator. We took the hyperparameter values 
from the photometric GP regression and used them as 
initial fit values for the RV GP kernel hyperparameters, 
with the exception of the hyperparameters h and cr, as 
the variation in the photometry and RV related to these 
quantities is not in equivalent units. We also place a 
Gaussian prior on 0 based on its posterior distribution de¬ 
termined from the photometric GP regression. We then 
apply a uniform offset to each set of RV measurements 
to remove RV zero point differences. 


The likelihood of the GP regression, given certain ker¬ 
nel hyperparameters and a Doppler amplitude value, de¬ 
termines the quality of the model for those given pa¬ 
rameters. This likelihood can be calculated as given in 
Equation 2, where the residuals r can be calculated as 


r = V 


i^sinl 


/ 27 r(t - tc) \ 
^ Pnrh ^ 


(4) 


where t is the vector of times of all measurements, v the 
vector of RV measurements, tc a time of transit, Porb the 
orbital period, and K the planetary Doppler amplitude. 
Sanchis-Ojeda et al. (2013) measured p to 10“^ of a day 
and and Porb to 10^ ^ of a day precision, so we can safely 
treat them as constants for this analysis. 


The best-fit GP kernel hyperparameters and the Ke- 
plerian Doppler amplitude are found via MGMG explo¬ 
ration of parameter space. We simultaneously fit the 
Keplerian planetary signal and a GP regression to each 
RV dataset with the Doppler amplitude and kernel hy¬ 
perparameters as free parameters in the MGMG chain. 
The emcee package contains an Afhne-invariant Monte 
Garlo Markov Chain Ensemble sampler, which evaluates 
the likelihood of the GP kernel hyperparameters to the 
measurement residuals a fter a Keplerian planetary sig - 
nal has been subtracted (Eoreman-Mackey et al. 2013). 
This is done for a plethora of steps in free parameter 
space. We draw the best-fit GP kernel hyperparameters 
and planetary Doppler amplitude as well as their errors 
from the posterior distributions generated through this 
MGMG exploration of parameter space, with 1-cr error 
corresponding to 68% confidence intervals in the MGMG 
posterior distributions. We plot the data and best-fit GP 
regression -1- Keplerian models in Eigurej^ 

After the MGMG posterior distributions have been cre¬ 
ated, we calculate a Gelman-Rubin statistic for each pa¬ 
rameter to ensure that the parameter chains have con¬ 
verged. Converg ence is deemed adequate when the G-R 
statistic < 1.01 (Gelman & Rubin 1992). In addition, 
we track the average and median log likelihood at each 
step in the chain. At first, the chains move to higher and 
higher likelihood parameter space. Once the best-fit solu¬ 
tion is found, the chains begin to take more and more un¬ 
likely steps, exploring parameter space more completely, 
and pushing the average likelihood below the median 
value. The step at which the average log likelihood be¬ 
gins dipping below the median log likelihood indicates 
the transition from burn-in period to the exploration of 
parameter space around the maximum likelihood solu¬ 
tion, and we remove the steps taken during burn-in from 


our analysis (e.g. Knutson et al.|2008). To check that the 
resultant fit is a good description of the data, we com¬ 
pute a x^ed best-fit model of each kernel variant. 

Einding X^^^ ~ 1 confirms that the quasi-periodic GP 
model and errors describe the data well. 


We allow as many of the MGMG parameters as possible 
be minimally constrained. We find that whenever the pe¬ 
riodic or quasi-periodic kernel is used, and the photomet¬ 
ric 0 value is not applied as a prior, the RV noise model 
is able to recover periodicity consistent with the stellar 
rotation period. This supports our assumption that the 
stellar surface features are the predominant signal seen in 
the RV measurements. Applying a Gaussian prior to this 
value allows us to better constrain all other MGMG pa¬ 
rameters. Physical boundary conditions were also placed 
on some hyperparameters to prevent the MGMG chain 
from moving into unphysical parameter space (see Table 
[^for details). 


4. RESULTS 


We find a Doppler amplitude of Kepler-78b of 1.86 ± 
0.25 m s“^, corresponding to a mass of 1.87lo.26 
using the quasi-periodic kernel with common temporal 
hyperparameters, excepting amplitudes and white noise 
terms, for the two RV datasets. We find that the quasi- 
periodic kernel provides the best fit to the Kepler pho¬ 
tometry based on visual signal inspection, and the valid¬ 
ity of the fit is confirmed through Xred calculation. Based 
on a Bayesian analysis, we find the strongest evidence for 
a quasi-periodic kernel regression of the radial velocity 
data over simpler kernels as well, and that stronger evi¬ 
dence exists for one set of temporal hyperparameters to 
describe both radial velocity datasets as opposed to two, 
suggesting that the dominant signal in both RV datasets 
is indeed astrophysica l. This supports the findings of 
Haywood et al. (2014). We provide the adopted model 
parameters for the best-fit quasi-periodic RV GP regres- 
3n -1- Keplerian orbit models chosen in Table 


4.1. GP kernel seleetion 

We find a planetary mass within l-a of our adopted 
mass result for all models tested. We make the assump¬ 
tion that the predominant signal observed in the photo¬ 
metric and both RV datasets comes from the star, and 
thus will modulate at the stellar rotation period. We 
find that when we model the photometric data with a 
periodic or quasi-periodic GP kernel, we recover a pe¬ 
riod hyperparameter consistent with previous estimates 
of the stellar rotation period, supporting our assumption 
about the predominant photometric signal. Similarly, 
when we model the RV datasets with a quasi-periodic 
GP with no priors, we recover a period hyperparameter 
consistent with the stellar rotation period or its alias. 
Thus, we conclude that the quasi-periodic GP regression 
is both well motivated and effective for describing the 
predominant signal in our datasets. 

We calculated BIG and AIC values for all of our RV 
models. We test the three kernels of interest both with 
shared hyper parameters as well as with fully indepen¬ 
dent hyperparameters in order to ensure that the varia- 
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Fig. 2. — RV of Kepler-78 versus time, measured by Keck-HIRES and HARPS-N. The HIRES and HARPS-N data are plotted as 
blue and red points respectively, with errors in RV shown. A Keplerian orbit signal with the calculated best-fit Doppler amplitude has 
been subtracted from the data. The colored lines correspond to the Gaussian process regressions with best-fit kernel parameters of the 
correspondingly colored RV measurements, where the shaded regions correspond to l-cr and 2-a uncertainties. Both datasets have been 
fitted with a single quasi-periodic kernel operator with common period, roughness, and lengthscale (0, ic, and A) hyperparameters, but 
separate covariance amplitude and white noise (h and a) parameters. 



Fig. 3.— The residuals of the HIRES and HARPS-N data after 
the quasi-periodic GP regression to the data with the Keplerian 
signal subtracted is removed, phase-folded at the known orbital 
period of the planet. HIRES data is shown in blue, and HARPS-N 
data is shown in red. The planetary signal model is shown by the 
dotted black line, and binned data are shown by black points. 


tion we observe is consistent between datasets and thus 
astrophysical, rather than local, in origin. Furthermore, 
we test our squared exponential kernel models with and 
without a white noise term for each dataset. We report 
these BIG and AIC values in Table [H 

We find that the quasi-periodic model we choose has 
AIG and BIG values with ABIG > —0.5 and AAIG > 8 
between it and all other models tested. A ABIG value 
of less than 2 suggests two models are indistinguishably 
likely, whereas a ABIG greater than 10 suggests strong 
evidence for the model with a lower BIG value. Interest¬ 
ingly, the BIG for the quasiperiodic and squared expo¬ 
nential kernels are almost identical, while there is signif¬ 
icantly less evidence for the periodic model (BIG of peri¬ 
odic model 18 points higher than SE and quasi-periodic 
models). Furthermore, the lowest BIG of all was found 
for the squared exponential kernel with two independent 
GP regressions, with a ABIC = —0.5 as compared to the 
chosen model. However, the AIC of the quasi-periodie, 
eoupled GP model is signifieantly lower than that of ei¬ 


ther of the other models, with a AAIC=8 between it and 
the next lowest AIC value, that for the quasi-periodie, 
uneoupled model, and AAIC =14 for the eoupled model 
with a squared exponential kernel. The relative likelihood 
of one model to another is given by e-^Aic/ 2 ^ indieating 
that the uneoupled quasi-periodie model is 1% as likely 
as the eoupled model, and the eoupled squared exponen¬ 
tial kernel model is 0.09% as likely. The independent 
GP, squared exponential kernel had a AAIC=19 when 
compared to the chosen model. We also tested removing 
the white noise parameter from the squared exponential 
model, finding that doing so raised the BIG by over 20 
points and the AIC by over 15. Thus, we determined 
that the white noise parameter was justified for all ker¬ 
nels. Strictly speaking, these information criteria require 
the assumption that the noise in the data is independent 
and identically distributed, which we know is not the 
case. However, since we use the same data for all BIG 
and AIC comparisons, they provide a valuable compar¬ 
ison between the different GP kernels. We report all 
relevant AIC and BIG values in Table O 


If we are modeling the stellar activity in both datasets, 
we would expect the hyperparameters for our simulta¬ 
neous models to be the same, with the exception of the 
non-temporal amplitude hyperparameter h and the white 
noise term cr, which would be different but related be¬ 
tween the models due to the different appearance of the 
disk of the star in different wavelength regimes. We ex¬ 
plored the relation of the GP kernel amplitude hyperpa¬ 
rameters, whose relationship we evaluated by calculating 
a, the ratio of the HARPS-N amplitude to the HIRES 
amplitude. The convergence of a to a best-fit value sug¬ 
gests that there is a steady relation in the ranges of the 
RV signals observed by HIRES and HARPS-N. Addi¬ 
tionally, t he value of a is bro adly consistent with the 
results ofiDesort et al.| (2007), who illustrate that the 


expected radial velocity contribution from starspots on 
a K dwarf is larger at shorter wavelengths, resulting in 
a larger signal in HIRES than in HARPS-N, as HARPS- 
N has significant wavelength coverage in the 600-700nm 
regime that HIRES does not. We show the distribution 
of a along with the other parameters of our MCMC cal- 
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TABLE 3 

Planetary Properties 


TABLE 4 

Adopted Quasi-periodic RV Model Parameters 


Property 

Value 

Name 

Radius, Rpi 
Orbital period Porb 

Kepler-78b 

1.20 ± O.O9R0 

0.35500744 ± 0.00000006 days 

Doppler amplitude, K 

1.86 ± 0.25 m s"! 

Mass, Mpi 


Density, Ppi 

e.olN g cm-3 

Iron fraction 

32 ± 26% 


Note. — - The name, radius, a nd orbital period values in Tableare 
taken from|Howard et al.|(|2013||. All other values have been calculated 
for this work. 


culation in Figure Any variation due explicitly to the 
correlated stellar activity should thus cause our models 
to have a shape (period and phase) that is consistent be¬ 
tween the two datasets. We allow the reader to confirm 
this visually in Figure 

Distributions for each parameter relative to all other 
parameters are shown in Figure for the quasi-periodic 
GP + Keplerian models. We see no clear correlation 
between Doppler amplitude and any other parameter 
tested, indicating that the Doppler amplitude observed 
is a real signal and not a systematic error introduced 
during our parameter fitting. 


We place a Jeffreys prior on the period hyperparame¬ 
ter 0 to weight shorter periods more heavily when fitting 
the quasi-periodic and periodic GP regr ession kernel hy¬ 


perpa rameters to the photometric data (Haywood et al. 


2014). When we do this, we are able to recover the ro- 


tation period of the star as the period hyperparameter. 
We then place a Gaussian prior on this hyperparame¬ 
ter in the RV GP + Keplerian model constraining it to 
the value and errors measured in the photometric data. 
We also note that although the coherence scale parame¬ 
ter w measured in the photometric data is equal to that 
measured in the RV data within errors, we do not place 
a similar Gaussian prior on this parameter, because we 
observe that it is correlated with both the amplitude and 
characteristic timescale hyperparameters, which were not 
necessarily related between the photometric and RV GP 
regression hyperparameter fits. This correlation between 
h, ^2, w and A visible in Figure is hard to interpret 
physically. We speculate that the positive correlation be¬ 
tween w and the h parameters may arise from a connec¬ 
tion to starspot size: as starspots grow larger, h grows, 
any white noise present will become relatively less impor¬ 
tant, and thus the function will become smoother over¬ 
all, increasing the w parameter. Similarly, the character¬ 
istic timescale A may also increase at large amplitudes 
for the same reason-as the uncorrelated noise at small 
timescales becomes less important, the lengthscale might 
be weighted more heavily toward larger values. Larger 
spots also persist longer on the stellar surface, as spots 


decay on timescales proportional to their size (Bumba 


1963). However, the fact that w was consistent tor the 


photometric and RV datasets whereas A was smaller for 
the photometric datasets despite their large differences 
in amplitudes may suggest otherwise. In addition, the 
amplitude parameters h and ^2 are not directly related 


Name 

Prior 

Value 

HIRES Amplitude h 

h>0 

ii.etls in s-i 

HARPS-N Amplitude /i2 

h2 > 0 

5.6): i'3 m S-I 

Amplitude ratio a 

- 

2 q+0.8 


(13.26-6>)2 

13.12 toXt. days 

Period 0 

g 2(0.12)2 

Coherence scale w 

- 

n oQ+o.os 

Characteristic timescale A 

A > 0 

26.1 tJi o days 

HIRES white noise crjitter 

a > 0 

2.lt°i m s-i 

HARPS-N white noise crjitter 

a > 0 

1.itS'.s m s-i 

Doppler amplitude, K 

- 

1.86 ± 0.25 m s"! 


to the range of the RV shift observed, and thus difficult 
to interpret physically. Further tests of these parameters 
are necessary in order to fully characterize their relation¬ 
ships, such as placing new priors on the coherence scale 
and characteristic timescale parameters to ensure that 
they are not biased by noise. 


5. INTERPRETATION AND DISGUSSION 


With our measured mass of 1.87lo!26 ^ radius 

(H13), Kepler-78b has a bulk density 


of 6.0l g cm suggesting a rocky composition similar 


of 1.20 ± 0.09 Ra 
LJi gem- 

to Earth {p^ = 5.52 g cm~^). Using t he tw o-component, 
rock-iron models from Eortney et al. (2007), we estimate 
an iron fraction of 32 A 26%, consistent wi th an Earth-like 
com positi on (iron mass fraction of 0.319, McDonough & 


Sun 


1995). These simplified models consider Kepler-78b 


as an iron core surrounded by a rocky mantel, and ac¬ 
count for compression that is important for higher planet 
masses. We ignore the effect of an atmosphere on the ra¬ 
dius of Keple r-78b due to its equilibrium temperature of 
1500-3000 K dSanchis-Ojeda et al.||M3 ). 


We now speculate about the implications of these mea¬ 
surements to the composition of Kepler-78b. In Eigurej^ 
we illustrate that the best-fit composition of Kepler-78b 
is indistinguishable from the composition of Earth. How¬ 
ever, Kepler-78b’s range in iron mass fraction stretches 
from almost purely rock (a Moon-like composition) to 
60% iron (a Mercury-like composition). This is consis¬ 
tent with a range of rocky solar system planets but dis¬ 
tinguishable from other solid celestial bodies (such as 
M-type asteroids, which are almost 100% iron by mass, 
or comets, made up of rock and ice with no iron at all). 
Since its best-fit composition and mass are more simi¬ 
lar to Earth than any other solar system body, Kepler- 
78b might have had a formation process that was similar 
to Earth’s. The proximity of Kepler-78 to its host star 
(0.009 AU) makes it unlikely that it formed in its current 
location, and migration to its current orbit is likely. 


The improvement in mass determination of this model 
relative to H13 is shown in Eigure Despite the fact 
that the mass measured by this method is a 6.5-cr mea¬ 
surement compared to H13’s 4-cr and P13’s 6-cr detection, 
the new measurement of Kepler-78b’s density is as pre¬ 
cise as those obtained by either of the two competing 
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Fig. 4. — Parameter distributions from our chosen GP + Keplerian model of the RV data plotted against each other. The median value 
of each parameter (blue lines) and 1-cr ranges (dashed lines) are shown. The parameter a = /i/h2. 


teams in 2013, who estimated its density to be 
(H13) and 5.57l^;3 g cm“^ (P13). Relative to P13, the 
errors on the density in this study are 15% smaller, al¬ 
though only marginally smaller than that of H13. This is 
because even though the new mass measurement is some¬ 
what more precise, the error on the density is driven three 
times more strongly by the error on the radius than the 
mass. Higher cadence photometry of the planetary tran¬ 
sit would allow for ingress and egress of the system to be 
observed, breaking the degeneracy between the impact 
parameter and transit depth. Estimates of the planet’s 
density could also be improved with more accurate mea¬ 
surements of the stellar mass, as converting RV measure¬ 
ments into a planetary mass is directly dependent on the 


stellar mass. 


Hatzes (2014) tested several traditional models on the 
same RV data analyzed here to check the robustness of 
the planetary detection and verify whether the planet 
could be found without previous transit knowledge. A 
range of Doppler amplitudes of the planet from 1.31 to 
1.96 m s“^ are reported, which is broadly consistent with 
the result of this work as well as that of H13 and P13. 


Hatzes (2014) reports the planet can be identified with¬ 
out prior knowledge of the system using a modified ver¬ 
sion of the parametric method originally tested by P13, 
the floating chunk offset model, to analyze both datasets, 
and reports a mass of 1.31 ± 0.24 M0, inconsistent at the 
1-cr level with the mass calculated in both this study and 
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TABLE 5 

Model MCMC Fit Diagnostics and Results 


Kernel function 

K (m s ^) 

BIC value (2GPs)'^ 

AIC value (2GPs)'^ 

Comment 

Quasi-periodic (adopted) ^ 

1.86 ± 0.25 m s'l 

1065.3 (1079.6) 

1039.2 (1047.0) 

with white noise, common le, 0, A 
but different h, a parameters 

Periodic ^ 

1.82 ± 0.29 m 

1083.8 (1090.6) 

1067.5 (1061.2) 

with white noise, common le, 0, 
but different h, a parameters 

Squared exponential 

1.92 ± 0.27 m s"! 

1066.4 (1064.8) 

1053.4 (1058.4) 

with white noise, common A but 
but different h parameters 


^ For the periodic and quasi-periodic kernel models, the period hyperparameter 0 is constrained by a Gaussian prior. This prior is shaped by the 
posterior distribution of the corresponding photometric GP hyperparameter. When this prior was not applied, the best-fit period found in the RV 
data was consistent with the stellar rotation period, but Doppler amplitude errors were larger. 

^ We measured the AIG and BIG values for all models with common RV GP hyperparameters and with 2 independent hyperparameter sets, which 
we call 2GPs here. We also measured BIG and AIG of the 2GP squared exponential kernel without a white noise term but do not include them in 
this table. 



Fig. 5. — A graph of mass versus radius for rocky planets. The 
constraints on Kepler-78b are shown by the red point. The es¬ 
timation of Kepler-78b by H13 is shown as the light gray error 
bars. Composition curves ranging from a pure water (blue) to pure 
iron composition (black) have been plotted. Earth-like (67% rock, 
33% iron) and Mercury-like (67% iron, 33% rock) compositions 
are denoted by green and brown curves, respectively. Solar sys¬ 
tem planets are shown as green squares. Other well-characterized 
exoplanets are plotted as black points. Exoplanet masses, radii, 
and their associated errors are from the Exoplanet Orbit Database 
(http://exoplanets.org downloaded on 23 August 2014). Planets 
with fractional mass uncertainties of over 50% are not shown. 


the previous studi es. The s pread in the Doppler ampli¬ 
tudes reported by Hatzes (|2Q14) illustrates that slight 


differences in the choice or noise model has a signifi¬ 
cant effect on the planetary signal extracted, and that 
the error on the planetary mass measurement has likely 
been underestimated. Thus, by starting with the sim¬ 
plest possible nonparametric descriptions of the data, we 
find a description of the data that is minimally complex 
yet strongly supported by Bayesian information criteria. 
Furthermore, the structure of the best fit GP kernel as 
well as the success of sharing kernel parameters to de¬ 
scribe both datasets supports our belief that we have 
robustly measured the astrophysical variability of the 


Kepler-78 system, and can obtain a mass measurement 
for Kepler-78b that is consistent with previous work. 

6. SUMMARY 

We performed a combined analysis of the photometric 
and RV data of Kepler-78 in order to better extract the 
RV signal of the planet Kepler-78b. We fit the data using 
Gaussian process regression, and test three simple kernel 
configurations. After testing multiple models, we find the 
strongest evidence for squared exponential kernel mod¬ 
els, and a quasi-periodic kernel model in which the hy¬ 
perparameters for the GP regression for each RV dataset 
are coupled (aside from the amplitude hyper parameter 
and the white noise term). We prefer the quasi-periodic 
model due to evidence indicated by AIG calculation and 
because it supports our understanding of the physical 
origins of the observed signals. We measure the Doppler 
amplitude and therefore the mass of Kepler-78b to 6.5- 
CF significance, comparable to or better than all previous 
mass measurements of this planet. We constrain the iron 
mass fraction of the planet to 32 ± 26%, illustrating that 
Kepler-78b is most likely Earth-like in composition. 

The analysis done in this work (and previous studies of 
this system) is possible because the orbital period of the 
planet Porb is an order of magnitude smaller than and 
not a harmonic of the rotation period Prot of Kepler-78. 
This made separation of the signals related to the stellar 
rotation and the signal related to the planetary orbital 
period possible. If the planetary period was a larger frac¬ 
tion or a harmonic of the stellar rotation period, decon¬ 
volving the signal due to the star and the signal due to 
the planet would be much more difficult. 

The true benefit of this analysis teehnique eomes from 
its nonparametrie nature. This analysis is particularly 
powerful in that even if the actual nature of the noise 
being modeled is unclear, a GP model can still be used 
to explore the nature of the noise and identify the most 
evident components. Due to the success of extracting a 
quasi-periodic signal at the rotation period of the star 
in independent datasets, we conclude that the predom¬ 
inant noise signal in the RV datasets likely comes from 
the stellar activity of Kepler-78. The variability of this 
noise, due to the growth and decay of active regions on 
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the stellar surface, cannot be fully parametrized with the 
information we have. Thus, testing the nonparametric 
Gaussian process noise model is a valuable exploration 
after earlier parametric and nonparametric models of the 
RV activity. Furthermore, the use of the Kepler photom¬ 
etry as a prior on our estimate of the RV activity makes 
our Gaussian process analysis especially valuable when 
photometric data is more readily available for a system 
than spectroscopic data, as is often the case. 

We plan to further test the robustness of this tech¬ 
nique by analyzing other RV datasets of exoplanetary 
systems. This method is particularly useful for analyz¬ 
ing stellar RV datasets over long time baselines which 
cannot be modeled easily because of the spot evolution 
on a timescale less than an order of magnitude larger 
than the stellar rotation period (as determined by the 
autocorrelation function analysis done by P13). Such 
scenarios can be easily described with a quasi-periodic 
GP. In addition, since kernel functions can be combined, 
any sort of physical combination of periodic and linear or 


exponential signals can be modeled with the GP, indicat¬ 
ing that it could be particularly useful in describing other 
noise modes seen in radial velocity studies as well as the 
already well-understood stellar surface signals. We hope 
to use this technique on a system with contemporaneous 
photometry and spectroscopy to explore the relationship 
between photometric and spectroscopic signals of exo¬ 
planetary systems. 
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